One-Stop Shopping for MBARI Upper Water Column Data : Analyzing signals from MARS using SSDS and Matlab

Following is an example of using Matlab and the loaddap tool to analyze high frequency data from the MARS2008 CTD instrument.

The code below will produce power spectra density plots from which signals may be discovered. We may also use this exercise to help us adjust sampling frequencies so proper experiment design and efficient use of MBARI's limited infrastructure.

% Load Original netCDF file into Matlab workspace.  URL may be discovered from the netcdf links on the web page:
% http://dods.mbari.org/data/ssdsdata/deployments/mars2008/current_qcPlots.html

loaddap('http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/mars2008/200706/mars2008_ctd1613_20081110_original.nc')

% May now do a 'who' to see what structures have been loaded.  Extract temp and time vectors from the SSDS data set:
temp = squeeze(Temperature.Temperature);
time = squeeze(Temperature.esecs);

% Construct regular 10-second interval time series for generating the spectral analysis
dT = 10;
time10 = [time(1):dT:time(end)];
temp10 = interp1(time,temp,time10);

% Apply 10% Cosine window to the time series and compute fft from demeaned temp time series
nt = length(time10);
rampUp = sin(linspace(0,pi/2,0.05*nt));
rampDown = cos(linspace(0,pi/2,0.05*nt));
nMiddle = nt - length(rampUp) - length(rampDown);
cos10 = [rampUp ones(1,nMiddle) rampDown];
temp_fft = fft(cos10 .* (temp10/mean(temp10)-1));

% Construct frequency axis and plot
f = linspace(0,1/dT,nt);
loglog(f, abs(temp_fft));
title('Time Series analysis of MARS CTD Instrument')
xlabel('Frequency (Hz)')
ylabel('PSD (Temperature)')
print -dpng temp_pds.png

Figure1. Power spectral density for Temperature at the MARS observatory. The peak in the power spectral density at 2.23 E-5 Hz is the M2 tide at 12.4 hours. Power density increases at high frequencies (periods less than 20 seconds) perhaps due to instrument noise.

pres = squeeze(Pressure.Pressure);
pres10 = interp1(time,pres,time10);
pres_fft = fft(cos10 .* (pres10/mean(pres10)-1));
loglog(f, abs(pres_fft));
title('Time Series analysis of MARS CTD Instrument')
xlabel('Frequency (Hz)')
ylabel('PSD (Pressure)')
 
% Use "peak<n> = ginput" in Matlab session to identify the peaks and label them with:
text(peak1(1), peak1(2) , sprintf('%.1f hours', 1/peak1(1)/3600), 'fontsize', 5)
text(peak2(1), peak2(2) , sprintf('%.1f hours', 1/peak2(1)/3600), 'fontsize', 5)
text(peak3(1), peak3(2) , sprintf('%.1f hours', 1/peak3(1)/3600), 'fontsize', 5)
print -dpng pres_pds.png

Figure 2. Power spectral density for Pressure at the MARS observatory. The sharp peaks at 25.8, 24.0, and 12.5 are due to the O1 (Principal diurnal lunar), P1 (Principal diurnal solar), and M2 (Principal semidiurnal lunar) tides. There is a shoulder of higher power density in the 10 to 1000 second band that we do not see in the temperature spectra.